import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns
import rasterio as rio
from rasterio.plot import show
from rasterio.plot import show_hist
import os
import fiona
!pwd
/c/Users/Bharath/anaconda3/envs/SG_Land_Reclamation_analysis
!ls
Approx_SG.kml debug.log Displaying images.ipynb Extracting_SG_from_Landsat8_images.ipynb Generated DataFrames HDB data Landsat_Image_processing_final_experiments.ipynb Landsat8_image_processing_initial_experimentation.ipynb Population_data_World_Bank.ipynb Population_density_and_Housing_density_analysis.ipynb Predicting the future of SG.ipynb SG_cropped_Landsat8_images SG_Landsat8_images SG_shapefiles
os.chdir('SG_cropped_Landsat8_images/')
!pwd
/c/Users/Bharath/anaconda3/envs/SG_Land_Reclamation_analysis/SG_cropped_Landsat8_images
!ls
2013 2014 2015 2016 2017 2018 2019 2020
os.chdir('2013')
SG_Bands_2013 = {}
for i in range(1,12):
var = 'B'+str(i)
img_file = 'SGimgB' + str(i) + '.TIF'
SG_Bands_2013[var] = rio.open(img_file, 'r')
os.chdir('..')
os.chdir('2014')
SG_Bands_2014 = {}
for i in range(1,12):
var = 'B'+str(i)
img_file = 'SGimgB' + str(i) + '.TIF'
SG_Bands_2014[var] = rio.open(img_file, 'r')
os.chdir('..')
os.chdir('2015')
SG_Bands_2015 = {}
for i in range(1,12):
var = 'B'+str(i)
img_file = 'SGimgB' + str(i) + '.TIF'
SG_Bands_2015[var] = rio.open(img_file, 'r')
os.chdir('..')
os.chdir('2016')
SG_Bands_2016 = {}
for i in range(1,12):
var = 'B'+str(i)
img_file = 'SGimgB' + str(i) + '.TIF'
SG_Bands_2016[var] = rio.open(img_file, 'r')
os.chdir('..')
os.chdir('2017')
SG_Bands_2017 = {}
for i in range(1,12):
var = 'B'+str(i)
img_file = 'SGimgB' + str(i) + '.TIF'
SG_Bands_2017[var] = rio.open(img_file, 'r')
os.chdir('..')
os.chdir('2018')
SG_Bands_2018 = {}
for i in range(1,12):
var = 'B'+str(i)
img_file = 'SGimgB' + str(i) + '.TIF'
SG_Bands_2018[var] = rio.open(img_file, 'r')
os.chdir('..')
os.chdir('2019')
SG_Bands_2019 = {}
for i in range(1,12):
var = 'B'+str(i)
img_file = 'SGimgB' + str(i) + '.TIF'
SG_Bands_2019[var] = rio.open(img_file, 'r')
os.chdir('..')
os.chdir('2020')
SG_Bands_2020 = {}
for i in range(1,12):
var = 'B'+str(i)
img_file = 'SGimgB' + str(i) + '.TIF'
SG_Bands_2020[var] = rio.open(img_file, 'r')
os.chdir('..')
def show_bands(dictionary):
Bands = dictionary
fig, ax = plt.subplots(3,3, figsize = (15,10))
num = 0
for i in range(3):
for j in range(3):
num += 1
key = 'B' + str(num)
ax[i,j].imshow(Bands[key].read(1))
ax[i,j].set_xlabel(key)
def experiment(dictionary, band_num):
Bands = dictionary
main_key = 'B' + str(band_num)
fig, ax = plt.subplots(3,3, figsize = (15,10))
num = 0
for i in range(3):
for j in range(3):
num += 1
key = 'B' + str(num)
if num == 8: # Excluding Band 8 since it has higher resolution and cannot be subtracted due to shape mismatch
continue
ax[i,j].imshow(Bands[main_key].read(1) - Bands[key].read(1))
ax[i,j].set_xlabel(main_key + " - " + key)
show_bands(SG_Bands_2013)
experiment(SG_Bands_2013, 5)
plt.figure(figsize = (15,10))
plt.imshow(SG_Bands_2013['B5'].read(1) - SG_Bands_2013['B4'].read(1) )
<matplotlib.image.AxesImage at 0x1f7221848b0>
show_bands(SG_Bands_2014)
experiment(SG_Bands_2014, 5)
plt.figure(figsize = (15,10))
plt.imshow(SG_Bands_2014['B5'].read(1) - SG_Bands_2014['B3'].read(1) )
<matplotlib.image.AxesImage at 0x1f7210e8160>
show_bands(SG_Bands_2015)
experiment(SG_Bands_2015, 5)
plt.figure(figsize = (15,10))
plt.imshow(SG_Bands_2015['B5'].read(1) - SG_Bands_2015['B3'].read(1) )
<matplotlib.image.AxesImage at 0x1f7221d2e50>
show_bands(SG_Bands_2016)
experiment(SG_Bands_2016, 5)
plt.figure(figsize = (15,10))
plt.imshow(SG_Bands_2016['B5'].read(1) - SG_Bands_2016['B3'].read(1) )
<matplotlib.image.AxesImage at 0x1f71f821d90>
show_bands(SG_Bands_2017)
experiment(SG_Bands_2017, 5)
plt.figure(figsize = (15,10))
plt.imshow(SG_Bands_2017['B5'].read(1) - SG_Bands_2017['B3'].read(1) )
<matplotlib.image.AxesImage at 0x1f71f70c490>
show_bands(SG_Bands_2018)
experiment(SG_Bands_2018, 5)
plt.figure(figsize = (15,10))
plt.imshow(SG_Bands_2018['B5'].read(1) - SG_Bands_2018['B2'].read(1) )
<matplotlib.image.AxesImage at 0x1f71d45d490>
show_bands(SG_Bands_2019)
experiment(SG_Bands_2019, 5)
plt.figure(figsize = (15,10))
plt.imshow(SG_Bands_2019['B5'].read(1) - SG_Bands_2019['B3'].read(1) )
<matplotlib.image.AxesImage at 0x1f71cddf0d0>
show_bands(SG_Bands_2020)
experiment(SG_Bands_2020, 5)
plt.figure(figsize = (15,10))
plt.imshow(SG_Bands_2020['B5'].read(1) - SG_Bands_2020['B3'].read(1) )
<matplotlib.image.AxesImage at 0x1f71f587910>
SG_imgs = {}
SG_imgs[2013] = SG_Bands_2013['B5'].read(1) - SG_Bands_2013['B3'].read(1)
SG_imgs[2014] = SG_Bands_2014['B5'].read(1) - SG_Bands_2014['B3'].read(1)
SG_imgs[2015] = SG_Bands_2015['B5'].read(1) - SG_Bands_2015['B3'].read(1)
SG_imgs[2016] = SG_Bands_2016['B5'].read(1) - SG_Bands_2016['B3'].read(1)
SG_imgs[2017] = SG_Bands_2017['B5'].read(1) - SG_Bands_2017['B3'].read(1)
SG_imgs[2018] = SG_Bands_2018['B5'].read(1) - SG_Bands_2018['B2'].read(1)
SG_imgs[2019] = SG_Bands_2019['B5'].read(1) - SG_Bands_2019['B3'].read(1)
SG_imgs[2020] = SG_Bands_2020['B5'].read(1) - SG_Bands_2020['B3'].read(1)
# Images that differentiate land and water the best
fig, ax = plt.subplots(len(SG_imgs)//2,2, figsize = (15,15))
num = 2013
for i in range(len(SG_imgs)//2):
for j in range(2):
ax[i, j].imshow(SG_imgs[num])
ax[i, j].set_xlabel(num)
num += 1
plt.tight_layout()
# Images flattened to 1D arrays
SG_imgs_flattened = {}
for i in range(2013,2021):
print(SG_imgs[i].shape)
SG_imgs_flattened[i] = SG_imgs[i].flatten()
print(SG_imgs_flattened[i])
(1216, 1877) [0 0 0 ... 0 0 0] (1216, 1877) [0 0 0 ... 0 0 0] (1216, 1877) [0 0 0 ... 0 0 0] (1216, 1877) [0 0 0 ... 0 0 0] (1216, 1877) [0 0 0 ... 0 0 0] (1216, 1877) [0 0 0 ... 0 0 0] (1216, 1877) [0 0 0 ... 0 0 0] (1216, 1877) [0 0 0 ... 0 0 0]
fig, ax = plt.subplots(len(SG_imgs)//2,2, figsize = (15,15))
num = 2013
for i in range(len(SG_imgs)//2):
for j in range(2):
ax[i, j].hist(SG_imgs_flattened[num], bins = 100, color = 'r')
ax[i, j].set_xlabel("Histogram for " + str(num))
num += 1
plt.tight_layout()
fig, ax = plt.subplots(len(SG_imgs),2, figsize = (15,15))
num = 2013
for i in range(len(SG_imgs)):
ax[i,0].imshow(SG_imgs[num])
ax[i,1].hist(SG_imgs_flattened[num], bins = 100, color = 'r')
ax[i,0].set_xlabel(num)
ax[i,1].set_xlabel("Histogram for " + str(num))
num += 1
plt.tight_layout()
#show_histogram(SG_Bands_2013)
# Finding the total number of pixels in the manually drawn shapefile
SG_Bands_2013_flat = SG_Bands_2013['B9'].read(1).flatten()
#type(SG_Bands_2013_flat)
total_pixels_in_cropped_img = (SG_Bands_2013_flat > 0).sum()
print(total_pixels_in_cropped_img)
1372139
# Thresholding
def Thresholding():
SG_binary_imgs_flat = {}
for i in range(2013, 2021):
SG_binary_imgs_flat[i] = np.zeros(len(SG_imgs_flattened[i]))
for j in range(len(SG_imgs_flattened[i])):
if SG_imgs_flattened[i][j] > 20000:
SG_binary_imgs_flat[i][j] = 1
return SG_binary_imgs_flat
SG_binary_imgs_flat = Thresholding()
SG_binary_imgs_flat_copy = SG_binary_imgs_flat.copy()
# Creating Binary image
SG_binary_imgs = {}
for i in range(2013, 2021):
SG_binary_imgs[i] = SG_binary_imgs_flat_copy[i].reshape((1216,1877))
fig, ax = plt.subplots(len(SG_imgs)//2,2, figsize = (15,15))
num = 2013
for i in range(len(SG_imgs)//2):
for j in range(2):
ax[i, j].imshow(SG_binary_imgs[num])
ax[i, j].set_xlabel("Binary Image: " + str(num))
num += 1
plt.tight_layout()
#Counting number of pixels of water
SG_water_pixels = {}
for i in range(2013, 2021):
SG_water_pixels[i] = (SG_binary_imgs_flat[i] == 1).sum()
#Counting number of pixels of land
SG_land_pixels = {}
for i in range(2013, 2021):
SG_land_pixels[i] = total_pixels_in_cropped_img - SG_water_pixels[i]
for i in range(2013, 2021):
print("Year: " + str(i))
print(" #Pixels: " + str(SG_land_pixels[i]))
Year: 2013
#Pixels: 774331
Year: 2014
#Pixels: 800777
Year: 2015
#Pixels: 791720
Year: 2016
#Pixels: 795768
Year: 2017
#Pixels: 801289
Year: 2018
#Pixels: 790573
Year: 2019
#Pixels: 798514
Year: 2020
#Pixels: 798597
# Dictionary that has the area of SG Land in squared km based on no. land pixels
SG_area = {}
# Data from https://data.gov.sg/
$C = constant$
$C * NumPixels = Area$
From 2015 data (Data from https://data.gov.sg/),
$C * 791720 = 719.1$
$ -> C = 0.0008947559682658218 $
C = 719.1/791720
for i in range(2013, 2021):
SG_area[i] = [C * SG_land_pixels[i]]
719.1/791720
0.0009082756530086394
print(SG_area)
{2013: [703.3059946698328], 2014: [727.3262525892992], 2015: [719.1], 2016: [722.776699843379], 2017: [727.7912897236397], 2018: [718.058207825999], 2019: [725.2708247865406], 2020: [725.3462116657404]}
!pwd
/c/Users/Bharath/anaconda3/envs/SG_Land_Reclamation_analysis/SG_cropped_Landsat8_images
os.chdir('..')
os.chdir('HDB data')
!ls
dwelling-units-under-hdbs-management-by-town-and-flat-type.csv estimated-percentage-of-singapore-resident-population-in-hdb-flats.csv estimated-resident-population-in-hdb-flats-by-town.csv Houses_by_type_2020.png metadata-estimated-resident-population-living-in-hdb-flats.txt metadata-number-of-residential-units-under-hdb-s-management.txt metadata-total-land-area-of-singapore.txt SG_population_data.csv total-land-area-of-singapore.csv World_bank_population_data.csv
SG_land_area_data = pd.read_csv('total-land-area-of-singapore.csv')
SG_land_area_data
| year | total_land_area | |
|---|---|---|
| 0 | 1960 | 581.5 |
| 1 | 1961 | 581.5 |
| 2 | 1962 | 581.5 |
| 3 | 1963 | 581.5 |
| 4 | 1964 | 581.5 |
| 5 | 1965 | 581.5 |
| 6 | 1966 | 581.5 |
| 7 | 1967 | 583.0 |
| 8 | 1968 | 584.3 |
| 9 | 1969 | 585.3 |
| 10 | 1970 | 586.4 |
| 11 | 1971 | 586.4 |
| 12 | 1972 | 586.4 |
| 13 | 1973 | 586.4 |
| 14 | 1974 | 587.6 |
| 15 | 1975 | 596.8 |
| 16 | 1976 | 602.0 |
| 17 | 1977 | 616.3 |
| 18 | 1978 | 616.3 |
| 19 | 1979 | 617.8 |
| 20 | 1980 | 617.8 |
| 21 | 1981 | 617.9 |
| 22 | 1982 | 618.1 |
| 23 | 1983 | 618.1 |
| 24 | 1984 | 620.2 |
| 25 | 1985 | 620.5 |
| 26 | 1986 | 621.7 |
| 27 | 1987 | 622.6 |
| 28 | 1988 | 625.6 |
| 29 | 1989 | 626.4 |
| 30 | 1990 | 633.0 |
| 31 | 1991 | 639.1 |
| 32 | 1992 | 641.0 |
| 33 | 1993 | 641.4 |
| 34 | 1994 | 646.1 |
| 35 | 1995 | 647.5 |
| 36 | 1996 | 647.5 |
| 37 | 1997 | 647.8 |
| 38 | 1998 | 648.1 |
| 39 | 1999 | 659.9 |
| 40 | 2000 | 682.7 |
| 41 | 2001 | 682.3 |
| 42 | 2002 | 687.1 |
| 43 | 2003 | 693.4 |
| 44 | 2004 | 696.2 |
| 45 | 2005 | 697.9 |
| 46 | 2006 | 699.5 |
| 47 | 2007 | 705.1 |
| 48 | 2008 | 710.2 |
| 49 | 2009 | 710.3 |
| 50 | 2010 | 712.4 |
| 51 | 2011 | 714.3 |
| 52 | 2012 | 715.8 |
SG_land_area = pd.DataFrame(SG_area, index = ["total_land_area"])
SG_land_area
| 2013 | 2014 | 2015 | 2016 | 2017 | 2018 | 2019 | 2020 | |
|---|---|---|---|---|---|---|---|---|
| total_land_area | 703.305995 | 727.326253 | 719.1 | 722.7767 | 727.79129 | 718.058208 | 725.270825 | 725.346212 |
SG_land_area = SG_land_area.transpose()
SG_land_area = SG_land_area.reset_index()
SG_land_area = SG_land_area.rename(columns ={"index": "year"} )
SG_land_area_data = SG_land_area_data.append(SG_land_area)
SG_land_area_data
| year | total_land_area | |
|---|---|---|
| 0 | 1960 | 581.500000 |
| 1 | 1961 | 581.500000 |
| 2 | 1962 | 581.500000 |
| 3 | 1963 | 581.500000 |
| 4 | 1964 | 581.500000 |
| ... | ... | ... |
| 3 | 2016 | 722.776700 |
| 4 | 2017 | 727.791290 |
| 5 | 2018 | 718.058208 |
| 6 | 2019 | 725.270825 |
| 7 | 2020 | 725.346212 |
61 rows × 2 columns
plt.rcParams.update({'font.size': 22})
plt.figure(figsize = (20,7))
plt.plot(SG_land_area_data['year'], SG_land_area_data['total_land_area'], '-o', label = 'Actual Area')
plt.plot(SG_land_area['year'], SG_land_area['total_land_area'], '-or', label = 'Estimated Area from Landsat 8 images')
plt.xlabel("Year")
plt.ylabel("Area in square km")
plt.legend()
plt.grid()
os.chdir('..')
os.chdir('Generated DataFrames/')
!pwd
/c/Users/Bharath/anaconda3/envs/SG_Land_Reclamation_analysis/Generated DataFrames
SG_land_area_data.to_csv('SG_land_area_data.csv', index= False)